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Abstract 

Recent developments in the theory of integrable models have provided the means of calculating 
dynamical correlation functions of some important observables in systems such as Heisenberg spin 
chains and one-dimensional atomic gases. This article explicitly describes how such calculations 
are generally implemented in the ABACUS C-|— |- library, emphasizing the universality in treatment 
of different cases coming as a consequence of unifying features within the Bethe Ansatz. 
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I. INTRODUCTION 

The Bethe Ansatz, introduced in 1931 [l| (very shortly after the appearance of quantum 
mechanics) is often viewed as an arcane subject, manifestly interesting for mathematically- 
minded theorists but not for a broad audience. At least two good reasons for this can 
De given. First, the subiect of integrability is rather specialized and demanding in itself 

iflnn 

21, 13|, U, 15| : it is after all made up of a collection of somewhat elaborate methods having only 
limited general applicability, and is therefore not typically considered worthy of inclusion in 
the standard condensed matter curriculum (the fact remains that the technology developed 
around integrable models has proven its usefulness in many ways, perhaps most importantly 
in providing many contributions to our understanding of strongly-correlated physics in one 
dimension 6j, but also in providing reliable beacons for the testing and benchmarking of 
more widely-used field-theory- or numerics-based methods for such systems). Second, and 
perhaps more importantly, is that the experimental relevance of integrable systems remains 
limited. Good realizations of one-dimensional integrable quantum systems do exit, but are 
the exception rather than the rule; when a good realization is found, integrability often 
cannot provide quantitative predictions for the most important observable quantities. 

This last difficulty is related to the long-standing inability of the Bethe Ansatz to provide 
results for even simple correlation functions of local physical observables. The source of this 
difficulty lies in the admittedly unwieldy form of Bethe wavefunctions, which are composed 
of a factorially large sum of free waves with nontrivial relative amplitudes. The action of 
local operators on these wavefunctions is difficult to write down, and since this forms the 
starting point in computing any correlation function, progress has been tedious, slow and 
limited. 

Recent years have however seen some important developments occur at a rapid pace in 
the calculation of dynamical correlations for finite systems at zero temperature. This stems 
from the appearance of a number of important results in the field of integrability. First and 
foremost is Slavnov's theorem on scalar products, providing initial results on form factors and 
correlation functions P, |8|]. The general problem of describing the action of local operators 
on Bethe wavefunctions was then resolved with the breakthrough solution of the quantum 
inverse problem for spin chains , opening the door to the computation of general 

correlation functions of the XXZ chain. The problem of obtaining correlation functions 
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from these representations was then initiated in ll|, ll2|, [iSj, with a study of traces over 



two-particle intermediate states. The ABACUS algorithm was thereafter developed in order 



to tackle the summation of all important general multipartic 
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e excitations in Heisenberg 



171]. The present article only 



spin chains and the one-dimensional Bose gas 
provides a 'conceptual-level' description of the features and functioning of ABACUS, only 
emphasizing unifying features of this class of calculations without going deeply into specific 
examples. My main objective in this article is to give a more or less detailed discussion of the 
strategies adopted in ABACUS, to define some terminology which will prove useful in future 



works, and to make some observations on the results obtained. A separate user guide [18 1 
will provide implementation-level details of the library and its usage, including important 
primary class and function declarations. Extensive and systematic results for various models 
will form the subject of specialized publications. The subject of integrability itself will be 
covered in a separate set of lecture notes [19i]. The general field of correlation functions of 
strongly-correlated one- dimensional quantum models is the subject of an upcoming review 



II. DEFINING THE PROBLEM 



Let us for definiteness consider a quantum-mechanical model on a lattice, with Hamil- 
tonian H which we assume for convenience (and without true loss of generality) to act in 
a Hilbert space formed by the tensoring of finitely many single-site local Hilbert spaces. 
Letting be the number of such sites and using j = 1, N as a site label, all operators 
can be reconstructed from the self-dual set Oj of local physical operators, where a is some 
on-site representation index. 

The general problem which we will address is the calculation of two-point zero tempera- 
ture equilibrium correlation functions of the form 

where (...) represents the ground-state expectation value and = {OjY, for arbitrary time 
differences t and distances j — j'. Equivalently, we can consider the Fourier transform of this 
correlation function, which (by a slight abuse of terminology) will be called the dynamical 
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structure factor (DSF): 




(2) 



A complete knowledge of the DSF for all momenta and frequencies naturally allows to 
reconstruct ([I]) without approximation. In reality, however, ([2]) itself is of more direct 
usefulness since it is often the quantity which can be related to energy- and momentum- 
resolved experimental measurements in linear response. 

Calculating the DSF for interacting models (by which is meant models where multiparti- 
cle states are not created by simple products of single-particle creation operators) is a very 
complex task. The product of quantum operators acting at different times and places is best 
dealt with by introducing a summation over intermediate eigenstates (labeled by greek in- 
dices), allowing to resolve the time dependence explicitly. Using the space Fourier transform 
O^. = Y2j G'^^'^Oj, summing over lattice sites and performing the time integration leads to 
the Lehmann series representation 



where Eq is the energy of the ground state |A°). The matrix elements {X\0\^) of local oper- 
ators in the eigenstates basis are also known as form factors (for continuum models); their 
norm square is also often called transition rate in the literature, in view of the correspondence 
of ([3]) to Fermi's golden rule. 

We can make a simple 'wish list' of the elements needed to provide a computation of such 
dynamical structure factors: 

1. an orthonormal eigenstate basis; 

2. form factors (matrix elements) {X\0^\n) of operators in this basis; 

3. a way to perform the summation over intermediate states. 

In general, this is an immensely complicated problem, which can typically be performed 
for models which correspond or can be mapped to free particles, or which benefit from a 
conformal field theory description. Simplistic treatments of interacting systems in general 
do not provide us with any of the three needed elements. 
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In one dimension, however, the technology of integrabihty does provide us with some of 
the three elements above, for a number of interesting models. The Bethe Ansatz provides 
the eigenstates basis, and the Algebraic Bethe Ansatz (together with the solution of the 
quantum inverse problem) provides us with form factors. The summation over intermediate 
states remains however an open problem, for which partial analytical results are possible 
in some restricted circumstances, but which still demands more often than not a numerical 
solution. 

The ABACUS library has recently been developed in order to perform the calcula- 
tion of dynamical structure factors of the known integrable models for which eigenstates 
and form factors can be explicitly constructed. ABACUS (Algebraic Bethe Ansatz-bascd 
Computation of Universal Structure factors) provides functions for representing and con- 
structing eigenstates, calculating form factors of a number of physical operators of interest, 
but also implements a scanning algorithm over intermediate states allowing to construct the 
correlations explicitly. My main motivations in developing ABACUS were primarily to pro- 
vide quantitative predictions first for experiments in low-dimensional magnetism and later 
on cold atomic gases, but also to provide results allowing cross-checking of more generally 
applicable alternative theoretical (analytical or numerical) approaches to the dynamics of 
strongly-correlated systems. 

The method which is presented here, and the results that it provides, undoubtedly suffer 
from a number of disadvantages and shortcomings. Among those, we can mention that: 

1. it is restricted to one-dimensional, Bethe Ansatz integrable models; 

2. it is restricted to the simplest known models, excluding all 'nested' systems; 

3. it only applies to zero temperature DSFs; 

4. it can only treat finite systems; 

5. it does not provide an exact, closed-form analytical expression, but only a numerical 
result. 

Restriction 1 is insurmountable: the whole ABACUS edifice begins and ends with the Bethe 
Ansatz. At best, what will be achievable is an extension of the framework to provide a 
description of correlations in perturbed integrable models. Restriction 2 is probably tem- 
porary, and exists simply because it is not yet known whether economical expressions for 
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form factors of local operators in nested systems exist or can be realistically found. It is 
also questionable whether restriction 3 is temporary or not: on the one hand, Bethe Ansatz 
can deal with finite temperatures quite straightforwardly for equilibrium expectation values. 
On the other hand, the problem of finite temperature Gibbs traces opens up a Pandora's 
box of complications within the ABACUS logic, and it is not yet clear whether these can be 
overcome. Restriction 4 comes from the fact that the eigenstates dealt with must be given 
by finite numbers of rapidities, and that the Hilbert space must have a measure of finiteness 
so that summing over the important intermediate states can be done in a finite time. It 
is thus not possible for ABACUS to give results in the strict infinite-size limit. The size 
dependence of the results is a complex affair, discussed in a later section. Restriction 5 is 
the most severe. ABACUS is and will remain numerical; whether simple exact closed-form 
expressions even exist for general DSFs is itself more than debatable at present. There are 
very good mathematical reasons to believe that such expressions cannot be found, except in 
various simplifying circumstances and/or limits. 

For completeness, and in order to relieve the discouragement of the reader, let us now 
fist the advantages of ABACUS: 

1. it works for very important 'canonical' strongly-correlated models such as Heisenberg 
chains and the Lieb-Liniger model; 

2. its implementation is relatively universal to the set of integrable models; 

3. it provides extremely accurate results for large systems; 

4. its accuracy is to a large degree energy and momentum independent. 

Advantage 1 is not strictly an advantage, but is worth emphasizing anyway: the models 
treated are strongly-interacting systems which do have experimental applicability, and an 
important role as beacons of reference for refining and fine-tuning alternative approaches. 
Advantage 2 comes from the fact that in the field of integrability, if one knows how to 
perform a 'trick' in one model, one can often apply the same trick to most if not all other 
integrable models. ABACUS is built from the start with the intension of exploiting this 
'universality' within the Bethe Ansatz. Advantage 3 is the most prominent: the accuracy 
level can be assessed with sum rules, and saturation levels beyond 99% are achievable for 
systems with many hundreds of particles, even in the more difficult limits. Relaxing this 
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requirement ever so slightly, say to 90%, one can go to systems with thousands of parti- 
cles. In any case, the system sizes achievable are comfortably much, much beyond anything 
which will ever be realistically achievable using exact diagonalization. Moreover, the system 
sizes attainable are also high enough that finite-size effects are drastically reduced, except 
perhaps in the vicinity of excitation thresholds, where the correspondence to singularities 
which develop in the infinite-size limit seems difficult to obtain. The fact that large systems 
can be treated is a consequence of the amount of preliminary information which is provided 
by integrability in the first place, and which can be built directly into the algorithm. The 
fact that ABACUS can give extremely accurate results for finite systems can itself be seen 
as an advantage, when applications to nano/mesoscopic quantum systems are considered. 
Advantage 4 is an interesting one, since many other methods (such as bosonization and con- 
formal field theory) are only valid at low energies. The common lore is that only 'universal' 
low-energy features should be computed by theorists, since these are the only ones which 
are unambiguously observable in experiments. A more careful assessment however teaches 
that such low-energy limits are precisely not accessible experimentally, and that reliable data 
read-out has to be done at higher energy scales. A less misleading statement about 'univer- 
sality' is thus that the extrapolation via scaling of the low-energy, universal result (which is 
not itself consistently observable) to somewhat higher energy scales (which are) can provide 
a good fit. However, the limitations of 'universal' approaches quickly becomes severe: in 
the case of spin chains, for example, inelastic neutron scattering provides data over the full 
Brillouin zone, showing clearly the band curvature and interaction effects which cannot be 
quantified using bosonization. ABACUS shows its strength here by being able to quantify 
these effects: it can compute the DSF at any energy scale, giving very precise data for any 
region of the momentum/energy plane. Since the Bethe Ansatz used by ABACUS gives 
exact wavefunctions irrespective of their energy, the accuracy of the data becomes more or 
less energy independent (since more types of excited states can live at higher energies, the 
accuracy there becomes limited by the restricted number of higher excited states that can 
be constructed; there thus remains an energy dependence to the accuracy, but it is weak). 
This feature of ABACUS makes it, at least in the mind of the author, the most appropriate 
method to give quantitative predictions for experiments. 

Let us now get down to the specifics. The bulk of the paper is organized as follows, 
mirroring the three elements of the above 'wish list'. Section IIIII provides an extensive 
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discussion of the conventions and terminology used in the classification and representation 
of eigenstates within ABACUS. Section [IV] then discusses how individual eigenstates and 
their norm are obtained, and how form factors of local operators are computed. Section |V] 
describes how the Hilbert space of intermediate states is scanned in order to achieve optimal 
results for the DSF. Section IVT] provides some example results obtained from ABACUS, with 
a discussion of some of their specific features. Finally, section IVIII presents a discussion of 
some generic features of the approach, and the paper ends with conclusions and perspectives 
for the future. 



III. CLASSIFICATION AND REPRESENTATION OF EIGENSTATES 



This section aims at providing a representation of eigenstates of Bethe Ansatz integrable 
models which is tailor-made for the purposes of calculating correlation functions. Such a 
representation will rely only on very generic features of Bethe eigenstates, so the discussion 
will not mention any model in particular (except to highlight exceptional aspects of specific 
cases) . 

Our integrable model H will have a Fock space JF of dimensionality dim(jF), which can 
be finite (for e.g. finite lattice models) of infinite (for e.g. continuum models with no UV 
cutoff). Typically, we will be able to easily separate this Fock space into a tensor product 
of spaces J^m with fixed number M of 'particles', with M being the charge of one of the 



trivially conserved quantities 



A. String basis; configurations 

The Bethe Ansatz gives us explicit wavefunctions for eigenstates of if in a fixed-charge 
subspace JF^,/, parametrized in terms of M rapidities A^, a = 1,...,M. In all generality, 
the rapidities of a given Bethe Ansatz solvable model live in the complex plane, and this 
introduces immense complications in the attempt to prove that the basis of Bethe eigenstates 
is complete. Under broad and reasonable assumptions, however, the common lore is that this 



basis is complete and that states can be understood and classified in terms of collections 



of (possibly deformed) generalized 'string' patterns 39|]. For a given model, let us call Ng the 



total number of possible string patterns. We let the label j run through the values 1, A^^ 
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with increasing j representing strings of increasing charge. In what follows, we call j the 
string level. A given eigenstate will have specific numbers Mj of strings of level j. We call 
a set {Mj}, j = 1, ...,Ns a base. The total charge M = '^jljMj where Ij is the charge of 
a string of level j is simply called the charge of a base, and we sometimes write the charge 
of a base as a subscript {i.e. {Mj}m)- The total number of strings, M^t = ^ M is 

called the string charge of a base. 

Since two-particle scattering preserves strings, the space J-'m further subdivides into sub- 
spaces ^{Mj}M having fixed charge-M bases. In other words, each partitioning of M into a 
charge-M base {Mj}M generates an independent Hilbert space. 

The exponential form of the Bethe equations is a set of transcendental equations for 
the rapidities. When written in logarithmic form, these equations can be interpreted as 
a mapping between the set of rapidities {A^} and a set of quantum numbers {/^} with 
j = 1, Ng and « = !,..., Mj at each level, 

NeUK) + E E ^sUK, AJ) = 27r/^ (4) 



k=l a=l 



where O^-^ is a kinetic function and 61^^^ is a scattering phase shift kernel. Since these equa- 
tions make use of the string hypothesis (of perfectly undeformed strings whose deviations 
are assumed exponentially small in system size), thus reducing the number of unknown pa 



rameters from the charge to the string charge, it seems historically appropriate |23l. l24j. |25| 
to call them the Bethe-Gaudin-Takahashi (BGT) equations, reserving the name Bethe equa- 
tions to the original full set (for a number of rapidity parameters equal to the charge) of 
coupled equations without any string hypothesis made. Although the rapidities are all in- 
terdependent, the quantum numbers are independent of one another, modulo applying a 
simple generalized Pauli principle stating that no two quantum numbers at the same level 
are allowed to coincide, i.e. for a ^ (3 E [l,Mj], we must have 7^ Vj. The quantum 
numbers thus provide us with the appropriate labelling of the eigenstates. We thus write 
|{/}) or I {A}) for our eigenstates, with the understanding that {A} are obtained from solving 
the Bethe-Gaudin-Takahashi equations with the choice {/}. 

The dimensionality of a given J^{Mj}m subspace coincides with the number of allowable 
choices of quantum numbers, since these obey a generalized Pauli principle {4^. These can 
be counted from the observation that choosing a base fixes the sets of possible quantum 
numbers {/^} which can be made at each level. In other words, the choice of a base {Mj} 



unambiguously defines a set {Z^,}; j = l,---,A^s of left- and right-limiting quantum num- 
bers Ilo,- aiid Iio,+ each level: any eigenstate in this base must have quantum numbers 
individually obeying the limiting inequalities G [— J^, _, -^4,,+], a = 1, Mj, j = 1, Ng 



4l| . The dimensionality of this subspace is thus 




dim(^{M,}J = lll ■ ). (5) 

7=1 V / 

We call a set of quantum numbers {/^} fulfilling the generalized Pauli principle and the 
limiting inequalities a configuration. A unique quantum number configuration is thus as- 
sociated with each individual eigenstate via the Bethe-Gaudin-Takahashi equations. One 
subtlety is that not all configurations lead to sensible string structures. If a set of quantum 
numbers {/^} is parity-invariant, then its associated rapidities will also be parity-invariant. 
If quantum numbers and /^^ for different levels ji and j2 within such a set are then 
simultaneously zero, both string centers A^^ and A^^ vanish according to the BGT equations. 
If the string charges are such that Ij^ = Ij^ mod 2, then the string hypothesis would pre- 
dict coinciding rapidities modulo exponentially small deviations {e.g. an overlapping 3- and 
5-string on zero rapidity leads to 3 pairwise superpositions). Since the string hypothesis 
patently fails in such cases, we call these states inadmissible, the usual states being admissi- 
ble. Inadmissible states do give proper wavefunctions involving (strongly) deviated strings, 
but the original Bethe equations must be solved. ABACUS ignores inadmissible states in 
the current implementation, since the numbers of inadmissible states necessarily are sup- 
pressed in a fashion as compared to admissible ones, and in any case their contribution 
to correlation functions can in some cases be shown to vanish identically 26 1. 

For the calculation of zero-temperature correlation functions as performed by ABACUS, 
the most important eigenstates are those for which Mi < M and Mj ~ 0(1) for j > 1, in 
other words eigenstates for which only a handful of 'higher' strings are present, and most 
rapidities sit at level j = 1. 



1. Representation of configurations in ABACUS 

Let us now describe how a configuration is represented and labeled in ABACUS. Let us 
start by considering level j = 1, the level at which we will make use of the most complicated 
representation. The Mi quantum numbers of this level are to be chosen within the set of 
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quantum numbers — _, — /^^ _ + 1^^, /^^ + — 1, /^^ +. We represent the possible choices by 
open circles, as per a) in Figure [TJ 42] In the same figure, b) represents the lowest-energy 
configuration, while c) and d) represent higher-energy configurations. The fact that b) 
represents the lowest-energy configuration comes by design of the rapidity parametrization, 
and can be ensured at all levels, for all models. 

We often use the term 'particle' to denote a given occupied quantum number, representing 
a string in the eigenstate. In most cases, we call 'excitations' of a configuration the holes 
within the base Fermi set (this being defined as the lowest-energy configuration for a given 
base), and the particles outside of it. Thus, configurations b), c), d) in Figure [H each have 
5 particles, but example b) has excitations, while c) has 4, and d) has six (this is not a 
strict rule: a simple exception are the two-spinon states of Heisenberg chains in zero field, 
where the spinous are dispersive 'holes' 27|, and the 'particles' living out of the base Fermi 



set are fixed (only one quantum number slot is available); since the positions of the two 
holes completely specify the state, we then say that such states contain two excitations). In 
circumstances where we need to be more specific, we will distinguish between 'dispersive' 
and 'non- dispersive' excitations. The latter have only one possible choice for their quantum 
number, and therefore are not part of an excitation family; the former, on the other hand, 
have many quantum number possibilities, and in most cases become true excitation branches 
in the infinite-size limit. 

a. Sectors. Excitations formed by a number of particle-hole pairs can be further sepa- 
rated into different classes, depending on which sides of the base Fermi set they live on. We 
first define four sectors as illustrated in Figure O Sector contains the quantum numbers 
obeying I > 0, I < Ip = ^ (right-hand side of the base Fermi set, including zero). Sector 1 
contains those obeying / < 0, / > —Ip (left-hand side of base Fermi set). Sector 2 contains 
the quantum numbers for which I > Ip and / < while sector 3 contains those obeying 

/ < —Ip and / > — /i, The numbers of excitations e|, i = 0, 3 at level j = 1, combined 



and written out as e\e\e\e\ is called the level 1 type [43|. The fact that excitations come 
from particle- hole pairs means that at level j = 1, we have the identity e\ + e\ = e\ + e\. 
Despite this slight redundancy, we still label the type as described above, since it then allows 
to visualize the overall quantum number excitation pattern of a state with a quick glance. 
Thus, in Figure [H the configurations b, c and d are respectively of level 1 type 0000, 1111 
and 1212 y. 
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FIG. 1: Representation of quantum numbers at level j = 1- In a), the set of possible quantum 
numbers at this level for a given base is shown as open circles. In this case, we must fill in Mi = 5 
of the possibilities defined by the limiting quantum numbers /<^^_ = 5 and = 6. b) shows 

the lowest-energy configuration, c) shows a relatively low-energy excited configuration with two 
particlc-hole pairs (four excitations) at this level, with the particles and holes staying close to the 
Fermi boundary. Finally, d) gives an example of a high-energy configuration with three particle- 
hole pairs (six excitations), with holes going further in and particles further out away from the 
Fermi boundary. 
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FIG. 2: The four sectors of excitations at level j = I. Sectors and 1 are respectively the right- 
and left-hand sides of the base Fermi set, while 2 and 3 are the right- and left-outsides of the base 
Fermi set. The excitation type for the upper configuration is 0000. The excitation type for the 
lower configuration is 1212 (see main text). 
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FIG. 3: The two sectors for levels j > I- The level j type for this configuration is 13. 

Levels j > 1 use a slightly less complicated representation. As mentioned before, in 
practice ABACUS only needs to construct states for which Mj ~ 0(1) for j > 1, and we 
therefore use only two sectors for these levels. As illustrated in Figure [3l sector at level j 
contains allowable quantum numbers P > 0, while sector 1 contains P < 0. The level j type 
is defined similarly to the level 1 type, but is here given by the simple combination c^Cq. 

The whole set of level j types for j = 1, Ns defines the type of the eigenstate. 

b. Mapping to Young tableaux. Within a given sector at a given level, we map the 
configuration of occupied and unoccupied quantum numbers onto a Young tableau. This is 
done differently in each of the four sectors of level 1, and in the two sectors of each level 
J>1- 

It is simplest to illustrate the general idea using a level j > 1 . This is done in Figure HI 
where a type 22 configuration is considered. Sector of that level has 2 excitations; there 
are six allowable quantum numbers. The lowest energy configuration of this type would 
have the two quantum numbers 1/2 and 3/2 occupied. This would be mapped to an empty 
tableau. Instead, in the configuration shown, the rightmost excitation sits on 11/2, four 
slots higher than its quantum number 3/2 in the lowest configuration. This displacement 
is represented by putting 4 boxes in the first row of the tableau. Similarly, the second- 
rightmost excitation sits on 5/2, two slots higher than its original 1/2. We therefore put 
two boxes in the second row of the tableau. In sector 1, we have two excitations, and 
we use the same logic but inverted right-to- left, the first row of the tableau in sector 1 
thus represents the left-displacement of the excitation which sits left-most in the lowest 
configuration. It should be more or less immediately clear to the reader that the generalized 
Pauli principle means that the tableaux that are thus obtained obey the Young tableaux 
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FIG. 4: Mapping quantum number configurations in a sector to a Young tableau, here for the two 
sectors at a level j > 1. 

rules. The tableaux in a given sector can also be read in a dual manner: if the row lengths 
represent e.g. the right-displacements of a set of particles/holes, then the column lengths 
represent the left-displacements of the associated holes/particles. 

At level j — 1, we need four tableaux. The sector tableau at level 1 represents the 
inward left-displacements (towards the middle) of the holes on the right of the base Fermi 
set. The sector 1 tableau similarly represents the right-displacements of the holes on the left 
of the base Fermi set. Sectors 2 and 3 represent respectively the right- (lcft-)diplacements of 
the excitations to the right (left) of the base Fermi set, similarly to the logic used in levels 

i>i- 

For a given sector i at a given level j, the number of rows in the corresponding tableau 
is equal to the number of excitations in this sector. The number of columns ric is given by 
the maximal possible displacement of the 'highest' excitation, which is simply the number 
of available quantum numbers nj in this sector minus the number of excitations. There are 

f 4 \ 

thus = I j ^ "^"^ different allowable tableaux (in other words, configurations) in 

this sector at this level. 

We number Young tableaux in the following manner: the empty tableau is given id — 0. 
We add boxes on the first row until the row is full, increasing the id by one for each added 
box. When we reach the maximal length, the next tableau is the one with a single box on 
the second row (and thus also on the first row). We then add boxes again on the first row up 
to the maximal length. The next tableau is then the one with two boxes on the second row. 
We proceed like this until both the first and second rows are full. The next tableau is then 
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FIG. 5: Young tableaux and their identification numbers. Tableaux 0-9 form the set of all ten 
tableaux with (rir, Uc) = (2, 3). Tableaux 0-19 form the set of twenty tableaux with {nr,nc) = (3, 3). 

the one with one box on the third row (and thus also on the second and first ones). This is 
done until the full tableau is obtained. An example of this is given in Figure [5l where we 
list all the (n^, ric) = (3, 3) tableaux and their identification numbers. It is simple to write a 
recursive algorithm providing the mapping to/from Young tableaux from/to identification 
numbers. 

c. Eigenstate labelling in ABACUS. We now have all the needed elements to describe 
the labelling of eigenstates within ABACUS. An eigenstate is uniquely identified by a config- 
uration. A configuration is uniquely defined by three integers: its base_id, its offset type_id, 
and its offset_id. 

The base_id is a long long integer specifying the set {Mj}. In general, we have Mi ~ M 
1, with Mj ~ 0(1) for j > 1. We thus write 

base_id = [##],,^J##]m,^_^...[##]m3[##]mJ#####]m„ 

in which [4i^4f\h[j signifies two integers giving M,-, e.g. 03 or 11 or other two-digit non- 
negative number less than 100. At level 1, in order to permit the labelling of states in 
large systems, we allow for a five-digit specification. The numbers should be read from right 
to left, since leading O's are scrapped by file or standard output. A type_id of 1040200293 
therefore means that the base is given by Mi = 293, M2 = 2, M3 = 4, M4 = 1 and Mj = 0, 
J > 4. 

The type_id of the offset of a configuration is represented as an integer 
type_id = e^^e^'e^^-'e^'-\...eleleielelelelel 
15 



(N.B.: the numbers should again be read /rom right to left). This is also represented as a long 
long integer in ABACUS, since in practice the code never puts more than nine excitations 
in a sector (allowing the type_id to be interpreted one character at a time), and uses only 
the lowest set of levels (meaning that the limitations on the length of integers do not affect 
the actual calculations performed). 

Finally, the offset_id is defined as follows. Recall that ci] is the total number of tableaux 
of sector i at level j. The tableaux identification number tj for this sector at this level thus 
runs through the values 0, 1, c?^ — 1. Given a set of tableaux identification numbers {t^} 
for each sector at each level, we thus define 

Ns 3 i-1 j-1 3 



j=l i=0 i'=0 j'=li"=0 

as the single integer defining all tableaux identification numbers. In other words, we order 

the (s, /) sector/level pairs as 



and the offset_id is obtained by multiplying each tableau id by the product of the number 
of possible tableaux at each lower sector/level pair. 

We therefore have defined a mapping between three integers, and the whole set of quan- 
tum numbers of an eigenstate: 

(baseJd, typeJd, offsetJd) {/^} 

where the string content and excitation pattern are clearly {i.e. humanly) legible from the 
first and second integers respectively. These numbers are used in all functions used by 
ABACUS, and in the data files that it produces. 

IV. CONSTRUCTING INDIVIDUAL EIGENSATES AND CALCULATING 
THEIR FORM FACTORS 

A. Solving the Bethe equations 

As described above, the Bethe Ansatz provides an economical pathway for obtaining 
eigenfunctions by replacing the diagonalization of the Hamiltonian by the process of solving 




(0, 0), (1, 0), (2, 0), (3, 0), (0, 1), (1, 1), (0, 2), (1, 2), (0, iV, - 1), (1, N,-l), 
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some sets of nonlinear coupled equations for one eigenstate at a time. The overall difficulty 
of solving the Bethe equations varies greatly from one model to the other. Also, within one 
model, the solution to the Bethe equations are easy to find for some of the eigenstates, while 
for some other eigenstates this might represent a very difficult challenge. 

The simplest situation is when all solutions to all Bethe equations are to be found in 
terms of real, finite rapidities. This is the case for example for the one- dimensional Bose gas 
(Lieb-Liniger model) in the repulsive regime: the proof is simple, and relies on the convexity 
the Yang- Yang action associated to the Bethe equations (jl]) on the field of real numbers 



281, see also details in [3]). For Heisenberg spin chains however, the corresponding Yang- 



Yang action is not convex, and complex solutions exist. It is a well-known fact that there 
is no, and cannot be, a good numerical method for solving coupled nonlinear systems of 
equations, so we have to input as much preliminary knowledge as possible in the procedure. 
The string hypothesis conveniently provides such preliminary information. 

In order to solve the Bethe-Gaudin-Takahashi equations, ABACUS mixes different meth- 
ods. The simplest method, used in the initial stages of the solution process, consists in 
performing some simple iterations. Defining the set {A;^^.^^} as the set of rapidities at iter- 
ation number i, we can proceed as follows. First, the initial set is defined as the solution 
to the decoupled conditions NQ-^^^i^X-'^ ^q-j) = 27rJ^. Subsequent sets are defined from the re- 
cursion relation Ar^L(^a,(i+i)) = ^^^^l - Y.k=i EJ=i ^scat(Ai,(i), (■)), possibly using some 
damping to help stabilize the process. Such simple iterations have the advantage that they 
often work, in the sense that some measure of convergence is achieved. Each iterative step 
is also rather quick, involving order operations. However, convergence is slow, especially 
when the solution is approached, and it is therefore desirable to accelerate the algorithm in 
various ways. 

One way to converge more quickly is to track the changes of the rapidities over a certain 
number of iterations. This gives a flow pattern to each rapidity: if this flow is sufficiently 
regular, it can then be extrapolated to an inflnite number of iterations in order to generate 
a new set of rapidities. In practice, this procedure of combining iterations (say four or flve) 
and extrapolations substantially accelerates convergence as compared to doing only simple 
iterations. 

A more elegant approach becomes available once the iterations have come sufficiently 
close to the true solution: the matrix Newton method. This has the advantage of converg- 
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ing quadratically, but has the disadvantages that 1) it is more computationally costly per 
iteration step, involving order operations, and 2) it is not always stable if the starting set 
of rapidities is too far from the actual solution (this instability is a possibility if the Yang- 
Yang action is not convex). Various run-time checks need to be made to see if this procedure 
is stable. One advantage of the matrix Newton method is that the Jacobian matrix which 
is needed, given by the matrix of derivatives of the left-hand side of dl]), is precisely the 
Gaudin matrix which is needed to compute the norm of the eigenstate (discussed in the 
next subsection), which also provides a slight acceleration of the algorithm. 

ABACUS considers the Bethe equations as solved provided a further iterative or Newton 
step gives a new set of rapidities for which the sum of square differences with the previous 
set is below a fixed convergence precision threshold, the latter being set to a fixed value 
somewhat larger than the available numerical accuracy (machine epsilon). In the case of 
purely real rapidities, this is the only convergence criterion applied. For states with strings, 
however, we have to be more careful. Namely, solutions to the Bethe-Gaudin-Takahashi 
equations do not always represent proper solutions to the corresponding Bethe equations, 
since the string hypothesis might not be verified to sufficident accuracy for that specific 
eigenstate. In spin chains, this occurs in many circumstances, including for eigenstates with 
strings having a very large rapidity, or when the magnetic field is near zero (for a more 
extensive discussion of deviated strings, see j^] and references therein). For eigenstates 
with string states, ABACUS therefore also computes the first-order string deviations by 
substituting the solution of the Bethe-Gaudin-Takahashi equations back into the original 
Bethe equations, yielding a new iteration for the (now full) set of rapidities, this time 
including their imaginary parts explicitly. If the sum of the string deviations (defined as 
the difference between the patterns obtained and pure, undeformed strings) is higher than a 
given fixed string precision threshold, the eigenstate is deemed unusable, and its contribution 
is excluded from the final result to prevent corruption of the data. If these deviations are 
small enough, the contribution from this state will however be kept. 

B. Calculating norms and form factors 

Once the set of rapidities {A^} of an eigenstate is known, it is then a matter of straightfor- 
ward number crunching to obtain the norm of the eigenstate. This is given by the Gaudin- 
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Korepin formula [29|, |30| in terms of the determinant of the Gaudin matrix, the latter being 
defined (as mentioned above) as the derivative matrix of the Bethe equations (in the case of 
string states, the norm can be similarly calculated from a reduced Gaudin matrix obtained 
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31|). The determinant itself is then 



from the Bethe-Gaudin-Takahashi equations, see 
calculated using LU decomposition. 

For form factors of local operators in non-nested Bethe Ansatz integrable models, the 
Algebraic Bethe Ansatz provides a representation in terms of a matrix determinant, similarly 
to the norm. In this case, however, the matrix depends on two sets of rapidities (one for 
each of the bra and ket states), and takes a different form depending on which operator 
is considered. The explicit formulas can be found in the literature. For the Bose gas, the 



density operator form factor is given in 0, Q] • The field operator form factor was given as a 
sum of determinants in |32j^, and simplified to a single determinant ir i |l7 |. For Heisenberg 
S =1/2 chains, the and operator form factors were given in jol. Iiol for the case where 
string states exist at most on one side of the bracket. The general expression, valid when 
either (or both) of the bra and ket states have strings, was given in [l5| 45|. Here again, 
the determinant is calculated by ABACUS using LU decomposition. 



V. SCANNING THE HILBERT SPACE 



Assuming that eigenstates can be classified, individually labeled and constructed as de- 
scribed in the previous sections, we now come to the more challenging third item on the 
introduction's 'wish list': the summation over intermediate states. 

For a given DSF in a given model, the perfect scanning algorithm will generate eigenstates 
in order of numerically decreasing absolute value of form factors. This is not trivial to achieve 
in the models we consider, since the intuitions we can develop relating to the solutions to 
the Bethe equations and the rapidity dependence of the form factors, remain incomplete. 

All is not lost, however. Four general but simple guiding principles can be used to define 
an efficient algorithm to scan the Hilbert space. The first principle is that i) for large 
enough system sizes, form factors obey an approximate continuity principle, in the sense 
that modifying an excited state in a small way {e.g. moving just one quantum number) 
does not dramatically change the form factor. This principle is ultimately related to the 
fact that form factors are analytic functions of the rapidities of the eigenstates involved, and 
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translates to the statement that if we find some excited state having a large contribution, 
we can exploit this by concentrating computational resources on the pool of states in its 
vicinity. 

The remaining three principles associate in turn to the three integers we use to label 
eigenstates. In general, we observe that form factors of local physical operators between 
the ground state and an excited state decrease in numerical value for increasing complexity 
of the excited state, by which we explicitly mean that form factors become smaller if ii) 
more/higher strings are included in the excited state, iii) more particle- hole excitations are 
created in the lowest level and iv) the excitations' quantum numbers are moved further and 
further away from the base Fermi set defining the lowest-energy state. 

These principles are not strictly true in all circumstances, and the ABACUS implemen- 
tation attempts to 'take them with a grain of salt' and probe sufficiently widely to gather 
all important contributions on the one hand, but not spend too much time on irrelevant 
contributions on the other. 



A. Climbing the Bethe tree 

The best way to visualize the scanning through the Hilbert space for the excited states is 



by using an analogy with another of Bethe's creations, namely a Bethe lattice [35|]. Starting 
from the lowest-energy state, one can then create particle-hole pairs in the lowest level, let 
these particles and holes move away from the Fermi configuration, add higher strings, let 
those disperse, and so on. 

Explicitly, the first excited state which is constructed is the lowest-energy state of the 
appropriate subspace, with the lowest base_id and zero type_id and offset_id. A particle- hole 
pair is then created, the typeJd becoming in turn 101, 110, 1010 and 1001. For each of these 
four branchings, a scan is made on the position of innermost (so sector or 1) excitation. 
If the contributions obtained to a chosen sum coming from these states is 'large enough' 
(whose meaning we will make explicit below), the next innermost excitation is 'raised' by 
a unit, and a new scan is made over the innermost excitation positions. The same logic 
is applied recursively in all sectors of all levels of the base in current use, as illustrated in 
figure El 

If a base, through this scheme, has given a large enough contribution, the base is then 
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FIG. 6: Climbing the Bethe tree for a given base. If an active excitation scan gives good enough 
results (in the sense defined in the text), a whole new generation of scans is initiated by descending 
the state (again, as explained in the text) in all possible ways, and starting a new active excitation 
scan for each such descendent. This procedure is continued recursively, and continues until good 
results are exhausted. Since each eigenstate is obtained via a unique such descendency path, the 
full Hilbert space can be faithfully probed. 

'complexified' at different levels, yielding a new family of bases on which the recursive 
scanning is also performed in a similar way until all large enough contributions from all 
bases have been exhausted. Very importantly, an ABACUS object called Scanned_lntervals 

maintains a map of the Hilbert space which shows which regions have been scanned, and what 
average results have been obtained per scanned region. This crucial run-time information 
will be used subsequently, as described in the next subsection. 

B. Optimizing the climb 

Defining the term 'large enough' used in the previous subsection leads us to the description 
of the highest-level procedure in ABACUS. The driving algorithm defines a real quantity 
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called running_threshold which defines a 'large enough' form factor as one whose value exceeds 
this threshold. The initial value of the threshold is chosen such that only relatively few (of 
order ~ A^^) intermediate states are constructed on a first run of the recursive Bethe tree 
climb. Once this first climb has terminated, a new lower value of running_threshold is chosen, 
and a new climb is initiated. 

The Scanned-lntervals object now fulfils its role, by 1) specifying to the scanning algorithm 
which parts of the Hilbert space of intermediate states can potentially yield contributions 
which are now 'large enough' according to the new running_threshold, and 2) which parts have 
already been scanned and should therefore be avoided to prevent state double-counting. 

When the second recursive Bethe tree climb terminates, the running_threshold is again 
reduced by an appropriate amount, and a third recursive climb is initiated. This is continued 
until a user-specified maximal allowed computing time is reached, before which ABACUS 
interprets and saves all the data to disk. Importantly, the Scanned_lntervals object is also 
saved, allowing a calculation which has successfully terminated to be restarted at a later 
stage, enabling to 'polish' results if required. A parallel implementation of ABACUS also 
exists, based on the same principles. This whole construction naturally provides an optimal 
use of the available computing resources. 

VI. RESULTS AND THEIR INTERPRETATION 

At the end of an ABACUS run, a large collection of form factor data becomes avail- 
able, composed of energy-momentum-form factor triplets, either in a file with additional 
state-by-state information, or conveniently binned into a form factor matrix (for large-scale 
computations producing immense numbers of entries, which could not all be individually 
saved to disk). The quality of this raw set of data (see an example in Figure [7j) can be 
measured by exploiting various sum rules, the simplest of which is to consider the sum- 
mation over all momenta and integration over all frequencies of equation [2l The resulting 
local expectation value {0°'0°') is in most cases known analytically. Since expression [3] is a 
summation over terms strictly greater than or equal to zero, the form factor data set can be 
summed up and compared to the analytical expectation value. A sum rule saturation per- 
centage can thus be assigned to each data set produced. Other sum rules can occasionally be 
used, for example the so-called /-sum rule which relates the integral over all frequencies of 
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Form factor results in a typical ABACUS run 




FIG. 7: Numerical value of (squared) form factors obtained during a typical ABACUS run (here, for 
the XX Z model at A = 0.6 and S^^^ = 0.1 for a small system of = 50 sites). The horizontal axis 
is simply the ordinality of the computed point. The continuous (in color: yellow) line represents 
a hypothetically-achieved perfect ordering of form factors in monotonically decreasing numerical 
value (an ideal algorithm would simply follow the yellow line). While still far from perfect (a factor 
of approximately 10 in speed could be gained by making the scanning perfect), ABACUS does 
follow such ideal lines reasonably closely in most circumstances. 

ujS"'°'{k, uj) to a known function of k. This is particularly useful since it allows to individually 
evaluate the quality of each separate momentum slice in the data set, which is beyond the 
reach of the integrated intensity. Sum rules are routinely checked by ABACUS: a summary 
file is produced, which contains information about the run. For each base and associated 
type, the number of scanned states is given, together with their contribution to the inte- 
grated intensity sum rule. Such summary files therefore provide a wealth of information 
about which types of excitations carry significant correlation weight, and which don't. 

Inevitably, for a given target sum rule saturation, the computation quickly becomes more 
intensive with increasing system size. The operation count for obtaining an individual state 
and its form factor scales like N^. In practice, the number of states that has to be included in 
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FIG. 8: Again for the typical example of the XX Z model at A = 0.6 and S^^^ = 0.1, the remainder 
of total integrated intensity sum rule after summing finite numbers of form factors is plotted. The 
number of states needed to achieve a required saturation increases more rapidly than polynomially 
in system size (as seen by the decreasing slope of the curves for increasing N), pointing to additional 
weak system size dependence in the ABACUS operation count exponent associated to the system 
size (see text for further discussion). 

the summation over intermediate states superficially seems to scale polynomially with as 
jY~3-5 (;iepending on the system and its parameters. As a rule of thumb, doubling the system 
size thus requires about two orders of magnitude more computational power to achieve the 
same saturation. In fact, the scaling with system size is a complex affair. If this scaling 
was of a well-defined degree, the curves of leftover sum rule weight as a function of number 
of contributions summed would have the same slope (assuming the efficiency of ABACUS 
to be system size independent). This is not the case, as can be seen from the example in 
Figure M and Table [H and there seems to be some additional (possibly logarithmic) system 
size dependence in the exponent. This is in fact quite natural: for larger and larger systems, 
we expect an increasing fraction of the correlation weight to be carried by states with higher 
excitation numbers. The system size scaling is discussed somewhat more in the discussion 
section, but a systematic study of it is beyond the scope of this article. 
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N = 50 


N = 100 


N = 200 


N = 400 


90% 


235 (4.9e-12) 


944 (6.8e-26) 


8772 (5.3e-54) 


117373 (3.5e-lll) 


95% 


385 (8.2e-12) 


2492 (1.8e-25) 


27096 (1.6e-53) 


573466 (1.7e-110) 


99% 


1515 (3.2e-ll) 


18490 (1.3e-24) 


469120 (2.8e-52) 




99.9% 


7478 (1.6e-10) 


141560 (l.Oe-23) 






99.99% 


26724 (5.7e-10) 


932706 (6.8e-23) 







TABLE I: Number of form factors needed to attain a given sum rule saturation level, for various 
values of system size. This is again for the example of the XXZ model at A = 0.6 and Sf^^ = 
0.1. The numbers in parentheses are the Hilbert space dimensionality fraction represented by the 
number of states mentioned, illustrating the efficiency of working in the chosen eigenstates basis. 
Missing numbers in the table could be obtained by simply running ABACUS for longer (all data 
here is from single-CPU runs). Although the number of states needed to achieve a fixed target sum 
rule saturation increases more rapidly than polynomially in as seen from Figure [71 the Hilbert 
space fraction needed also decreases. 

The finite system dynamical structure factor which ABACUS calculates, as is clear from 
equation [3], is not a smooth function in the energy-momentum plane, since it is simply com- 
posed of isolated delta peaks of various heights, aligned on well-defined lines in momentum 
but distributed in a very irregular pattern in energy (due to the adopted periodic boundary 
conditions on the wavef unctions, momentum is always an integer multiple of 27r/A^; energies 
are on the other hand certainly not equispaced in our interacting models). Only in the ther- 
modynamic limit does the DSF become a continuous function; to obtain this, some form of 
smoothing of the ABACUS results is necessary. Since the momentum is already sitting on 
a regular lattice, only the energy delta function is smoothened into a gaussian. This has 
the unfortunate effect of blurring some would-be sharp excitation thresholds, but allows to 
obtain easily interpretable density plots of the DSF. The width of the gaussian is chosen 
so as to be somewhat larger than the typical 2-excitation energy level spacing, so that the 
delta functions blend into a smooth function representing the density of states in this region 
of the Brillouin zone. Since this 2-excitation energy level spacing is typically of order 1/A^, 
the negative effects of the gaussian quickly become less significant for bigger systems. 
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VII. DISCUSSION 



A. Size dependence of contributions 

The dependence on system size of a given DSF in a given model is very complex. First 
of all, it may take different forms depending on which model parameters are used, and on 
which part of the energy-momentum plane one is studying. Changing system size affects i) 
the number of eigenstates in the Hilbert space, ii) the relative number of states of different 
bases, iii) the quality of numerical solutions to the Bethe equations, in particular the relative 
number of acceptable and unacceptable quality string states that can be constructed, iv) 
the form factor of a given identifiable state, and finally v) the positional distribution of 
eigenstates in energy (density of states). 

Point i) is trivial: the dimensionality of the Hilbert space is factorial in system size, 
and poses a severe limitation on the attainable size. Point ii) is more subtle and inter- 
esting. Namely, the number of n-excitation states scales approximately as (system size)"; 
therefore, the number of 2,3,4, ... excitation states scales differently, families with more ex- 
citations having an increasingly larger number of members. For any n-excitation state, we 
can generically find order (size)^ states with n + 2 excitations, and determining how these 
different bases relatively contribute to the final answer for large systems is thus an extremely 
complicated problem. 

We can however make some general statements. Suppose for definiteness that our ex- 
citations only come in pairs (as in e.g. the Lieb-Liniger model). For very small size (few 
particles), 2-excitation states will carry a fraction C2 ~ 1 of the sum rules. As the sys- 
tem size N increases, we will see the 2-excitation contribution 'leak' (as guaranteed by the 
preservation of sum rules) into the increasingly more numerous 4, 6, 8... -excitation ones. C2 
will be strictly monotonically decreasing in N, e.g. ^ < ^N. Despite this, the limit 
limAT^oo C2 might remain finite (as is the case for spinous in the zero-field isotropic or gapped 
Heisenberg antiferromagnet; for the gaplcss region, it seems reasonable for cj for j finite to 
vanish in the thermodynamic limit, except at zero field). C4, on the other hand, will not be 
monotonic: it will increase for small N but eventually reach a maximum at some inflexion 
value N4, after which it wiU decrease, i.e. ^ > 0, N < N4 hut ^ < 0, N > N4. The 
limiting value limjv^oo C4 might also remain flnite. There will similarly be increasingly large 
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inflexion values Nq, Ng, A^4 -C Nq ^ Ng...., but the determination of these inflexion 
values is beyond the reach of the current implementation. 

This scenario, although more or less inevitable in the mind of the author, remains con- 
jectural since even the lowest of these inflexion points (A^4) seems to lie above currently 
achievable system sizes. In any case, we can expect the finite-size result for large enough 
systems to closely resemble the infinite-size result (except perhaps in the immediate vicinity 
of excitation thresholds): while each term in the series C2 + C4 + C6 + ... might have substantial 
size dependence, the sum itself {i.e. the full DSF) is seen to have much smaller size depen- 
dence: increasing contributions (as system size increases) from higher excitation numbers 
tend more or less to compensate the decrease of the lower excitation number contributions 
(at least in gapless models), since these two types of contribution resemble each other. Thus, 



for example, the contribution of four-spinon states in the Heisenberg magnet |33| looks like 
a rescaled two-spinon contribution, and we can expect more or less the same of yet higher 
excitation numbers (this, disregarding the different high-energy tails, which carry very little 
correlation weight anyway). 

VIII. CONCLUSION AND PERSPECTIVES 

Much remains to be done in the general field of dynamical correlation functions of inte- 

□ 

grable models, and the reader is referred to [20| for an extensive discussion of this subject. 
As far as ABACUS is concerned, various improvements of the sub-algorithms for solving 
Bethe-Gaudin-Takahashi equations are possible, as well as further optimization of the state 
scanning algorithm. A better handling of deformed string eigenstates is also necessary: this, 
at the moment, represents the most severe limitation of the applicability of ABACUS for 
some DSFs of spin chains in small magnetic fields, but requires a rather more elaborate 
treatment to be properly dealt with. For certain for example gapped antiferromag- 

nets, the correct state counting itself is not fully known for arbitrary bases due to difficulties 



in defining the quantum number limits (see 3J| for further discussion of this point). 



The restriction to zero temperature would also need to be overcome, possibly by using 
some ideas from the thermodynamic Bethe Ansatz. On a more theoretical front, represen- 



tations of form factors for nested systems such as the Yang permutation model 361] or the 



Hubbard model l37| are needed before ABACUS can be used. Once these are available 
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however, the 'universal' ABACUS logic will be quickly implementable for these systems as 
well. 

d. Data availability Requests to the author for calculations of dynamical structure 
factors for specific models and values of the parameters are welcome. An online database of 
ABACUS computations will also eventually be made available. Alternately, the C++ source 
code can be obtained from the author on request. 
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